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Abstract 

Discrete choice models are commonly used by applied statisticians in numerous 
fields, such as marketing, economics, finance, and operations research. When agents 
in discrete choice models are assumed to have differing preferences, exact inference 
is often intractable. Markov chain Monte Carlo techniques make approximate infer- 
ence possible, but the computational cost is prohibitive on the large data sets now 
becoming routinely available. Variational methods provide a deterministic alternative 
for approximation of the posterior distribution. We derive variational procedures for 
empirical Bayes and fully Bayesian inference in the mixed multinomial logit model of 
discrete choice. The algorithms require only that we solve a sequence of unconstrained 
optimization problems, which are shown to be convex. Extensive simulations demon- 
strate that variational methods achieve accuracy competitive with Markov chain Monte 
Carlo, at a small fraction of the computational cost. Thus, variational methods permit 
inferences on data sets that otherwise could not be analyzed without bias-inducing 
modifications to the underlying model. 



1 Introduction 



Discrete choice models have a long history in statistical analysis, appearing in applications 



as varied as the analysis of consumer choice data ( |Guadagni and Little||1983| |Fader and 

Hardie] [T9"9~6 ), transportation planning ( |Theil 1969, McFadden| |1974 l Ben-Akiva and Ler- 

man||1985| ), economic demand estimation ( |Train et aL|1987| |Revelt and Train|fl998[ ), new 

product development ( |Moore et al.| [T999), portfolio analysis ( |Uhler and Cragg|[1971[ ) and 
*The authors gratefully acknowledge research assistance from Sandy Spicer and Liz Theurer. 



1 



health services deployment (Hal l et aLl|2002[ ). They apply to situations where agents (also 



called choosers or decision-makers) select items from a finite collection of alternatives (the 
choice set), either once or repeatedly over time. For example, in a marketing context, agents 
are "households"; each household makes a number of "trips" to a store, and we observe the 
items selected for purchase on each trip. 

Heterogeneous discrete choice models, which allow preferences to differ across agents, 
are based on a hierarchical regression formulation. We have agents numbered h = 1 , . . . , H, 
each with an unseen parameter vector Oh encoding preferences over item attributes. We ob- 
serve one or more choice events Yh ~ p(jh I Oh) per agent. The Oh's are modeled as 
independent draws from a prior distribution p(0h \ (ft), where <ft is a hyperparameter. This 
prior represents the heterogeneity of preferences across the population. Inference in such 
a hierarchical model allows us to pool information across decision-makers. If we use an 
empirical Bayes point estimate of eft (Robbins] |1955[ ), the posterior distribution of each $ h 



depends on all of Y\, . . . , Yh, through the common estimate (j>. In a fully Bayesian setup, 
integrating out the random variable (ft creates similar dependence. 

The marginal likelihood corresponding to one agent in a heterogeneous model is 

p(yn I 4>) = J p(yh I o h )p{O h I (ft) dO h . (i) 

In most cases, including the "random utility" discrete choice model we study in this pa- 
per, ([T|) does not exist in closed form. As a consequence, we must use approximate methods 
both for empirical Bayes and fully Bayesian inference. 

A standard empirical Bayes technique is to approximate ([T]) using Monte Carlo inte- 
gration. But to match the asymptotics of maximum likelihood, the number of draws per 
agent must grow faster than the square root of the number of agents (Train 2003[ ), which 



is infeasible for large-scale problems. The usual approach to the fully Bayesian random 
utility model is Markov chain Monte Carlo (MCMC) simulation ( |Albert and Chib|[T993 
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Allenby and Lenk||1994| |Allenby and Rossil|1999[). MCMC provides approximate draws 



from the joint posterior distribution on 0\, . . . , Oh and <p. The draws enable the estimation 
of ([T]) and related integrals. However, the more agents there are in the data set, the more 
MCMC output we need to collect and store - even if we are only interested in 0, we still 
need repeated draws for all of 0\ , . . . , Oh- 



Variational methods (Jordan et al. 1999; Wainwright and Jordan 2003) offer a deter 



ministic alternative for approximate inference. With variational inference, we maximize 
a data-dependent lower bound on the marginal likelihood ([T]), over a set of auxiliary pa- 
rameters distinct from the model parameters. In the fully Bayesian specification, the end 
result is an approximate joint posterior distribution for 0\, . . . , Oh and cp. For empirical 
Bayes, variational techniques lead to a point estimate 4> as well as an approximate posterior 
distribution for the #/ 7 's. 

The main advantage of variational methods versus MCMC is computational efficiency. 
Variational inference algorithms typically converge to their final approximation in far less 
time than it takes to generate an adequate number of MCMC draws. This advantage comes 
at the cost of a biased approximation, in contrast to the consistency guarantees that accom- 
pany MCMC. We give evidence in Section [4] that, for our random utility discrete choice 
model, variational bias is negligible, and the computational speedup is very large. Varia- 
tional convergence is also easy to assess, in contrast to MCMC. 

Furthermore, the size of the variational representation is fixed, while the size of the 
MCMC approximation increases with the number of draws. Variational techniques can 
therefore be applied to much larger data sets. For example, with 10,000 decision-makers 
and 20,000 total MCMC draws, using a 25-dimensional 6h and corresponding cfr, the MCMC 
representation of the posterior exceeds 2 GB — if we discard half the draws for burn-in of 
the chain. In fact, many data sets today contain observations on millions of agents, mod- 
els can contain far more than 25 preference parameters, and MCMC chains may require 



hundreds of thousands of iterations. 

These difficulties are well known. MCMC is rarely applied to large-scale heterogeneous 
models. Indeed, to address scalability, it is common to work with data from a subset of 
individuals, or a subset of choice items. However, this approach discards information that 



is valuable in the inferential process, and it can lead to biased estimates (Bradlow and 
Zanutto|2006l ). 



In this paper, we derive variational algorithms for a common discrete choice model - 
the mixed multinomial logit (MML) model. We study this model because it the workhorse 
of discrete choice theory and is well known in many disciplines, including economics and 
marketing. There are other popular discrete choice models in the literature, but the mixed 
multinomial logit has appeal: it is conceptually simple, yet still exhibits the MCMC infer- 
ence issues just described. 

The rest of the paper is organized as follows. Section [2] presents the details of the 
MML model. In Section[3j we describe variational procedures suitable for empirical Bayes 
and fully Bayesian inference under the MML model. These procedures include a novel 
application of the delta method for moments to variational inference. In Section |4} we 
compare variational methods to MCMC for the MML model, using an extensive suite of 
simulated data sets. Section [5] closes with discussion and future directions. Technical 
arguments and derivations are relegated to Appendixes |A} [Bj and[Cj 

2 The mixed multinomial logit model of discrete choice 

Let there be H agents, indexed h = 1, . . . , H. We observe a total of Th choice-event 
outcomes for agent h. At each choice event, the agent selects from among a fixed set of 
J items, indexed j = 1, . . . , J . The items are differentiated according to K attributes, 
indexed k = 1, . . . , K. The j'th item's value for the kth attribute can vary across agents 
and from one choice event to another. For example, households might shop at different 
stores charging various prices for the same good, and the price of a good may change over 



time within a single store. We denote by Xht the J x K matrix of attribute values, also 
called covariates, that agent h encounters at her fth choice event. The jth row of Xht is 
denoted xj t j . The outcome of this choice event is the observed categorical random variable 
jht, which we represent as a J x 1 indicator vector. 

We use the observed (x^, Jht) pairs to infer which attributes have the strongest as- 
sociation with item choice. To this end, let Uhtj denote the utility that accrues to agent 
h if she chooses item j at her fth choice event. This approach, called a "random util- 
ity model" ( Train||2003 ), assumes utility is a noisy linear function of the attributes: Uhtj = 



P~h x htj + eht j • Here, /?/, is a K x 1 vector of agent-specific "tastes" or "preference loadings" 
for the item attributes, and ehtj is a random error term representing unobserved utility. 

We assume that each agent, at each choice event, selects the item maximizing her utility. 
In the mixed multinomial logit model, we further assume the random error terms ehtj are 
iid from a Gumbel Type 2 distribution. The implied choice probabilities turn out to be 

Piyhtj = 1 I x ht , Ph) = ^ 7~nT 7 ' J = h...,J (2) 

2Lj' ex-Vifih x htj') 



( McFadden] 1 1 974[) . In discrete choice modeling, the right-hand side of ([2]) is called the 



"multinomial logit" distribution, denoted MNL(%,^/,). It is essentially the same as the 
multi-logistic function used in polychotomous logistic regression, and it is often called the 
soft-max function in machine learning research. 

We further assume that are iid from a ^-variate normal distribution with mean 
vector £" and co variance matrix Q, which we write as Nk(£, F° r empirical Bayes 
estimation, the model is now completely specified: 



Jht I x ht , Ph ~ MNL(x fa , p h ), h = l,...,H, t = l,...,T h , (3) 
ft h \C,n ~A^(c,Q), h = l,...,H. (4) 



5 



The top-level parameters f and Q, to be estimated by maximum marginal likelihood, rep- 
resent the distribution of attribute preferences across the population. In particular, Q. gives 
us information about the correlation of preferences between agents. 

A fully Bayesian approach requires hyperprior distributions for £" and Q.. As is standard, 
we use conditionally conjugate distributions: 



In ([5]), fa and Qq are pre-specified hyperparameters; W~ l (S~ l , v) is the inverse Wishart 
distribution with scale matrix 5 _1 and v degrees of freedom; and S and v are hyperpa- 
rameters fixed in advance. We call the fully Bayesian approach to MML model inference 
"hierarchical Bayes." 

3 Variational inference for the MML model 

We have presented the component hierarchical distributions in the mixed multinomial logit 
model. Now we turn to the question of estimation and inference procedures. In the follow- 
ing, variable names inside of /?(•) are used to distinguish among densities: we denote the 
pdfs in ©-© by p(y ht \ x ht , fi h ), p(fa \ £, Q), p(C \ fa, Qo)» and p(Q \ S, v), respec- 
tively. We let V = { Xht, yht} denote all observed variables, i.e., the data. 

For the empirical Bayes version of the MML model, the posterior density of the latent 
preference vectors, p{fa.H I f , Q), is 



C I fa, ~ Af K (j3o,£lo), 



Q\S,v 



(5) 



H 



p(p h \c, n)YlZ\p(yht\xht,Ph) 



n 



/ p(fi h | C, Utl PiJht I x ht , fa) dfa 



(6) 



The joint posterior density for hierarchical Bayes, pifa.H, <T> ^ I ^)» 1S 



p(o nf=i p(fa i c, n;=i po^ \ x ht , fa) 



(7) 



/ p(c) p(Q) ULi Jpifa l c, nil p(yht I x ht , fa) dfa, d{ da 
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The numerator in both cases is the joint density of latent and observed variables, computed 
by multiplying together the densities defined in the model hierarchy @-([5]). 

The integrals appearing in these posterior densities have no closed form. As a conse- 
quence, exact inference is intractable. Variational inference is a deterministic alternative to 
the MCMC methods usually applied to this problem (Ros si et al.|2005| ). A variational algo- 



rithm selects from a pre-specified family of distributions Q the best approximation to the 
true posterior distribution. We define Q so that all of its members permit tractable probabil- 
ity calculations. Then, wherever we need the true posterior, such as for an expectation, we 
use the approximating variational distribution instead. This plug-in idea underlies MCMC 
methods as well - in place of the true posterior, we substitute the empirical distribution of 
the MCMC posterior draws. 



3.1 Variational empirical Bayes 

In this section we give an overview of a variational algorithm for approximate empirical 
Bayes estimation in the MML model. Appendix |A~| fills in the details. We first specify 
a family of approximating distributions Q : = {qifti.H 1 1) : 1 e A) for the true pos- 
terior distribution (|6]). Since this posterior factors over h, we take Q to be a family of 
factored distributions as well, so that qifiuij I ^) : = Tlhlifih I ^h)- In particular, each 
factor q(fih I ^h) is a A'-variate normal density, with mean fih and covariance matrix E/,. 

For the particular data set at hand, we want to find qifti.H I A.*), the best approximation 
in Q to the posterior distribution p(P\-.H I £, Q.). To make the idea of a best approxi- 
mation precise, we measure discrepancy with the Kullback-Leibler (KL) divergence (also 
called the relative entropy). Shortening q(fii : H I ^) to qx, the optimal variational parame- 
ters are given by 

X* = argmin KL \qx \ \p\- (8) 
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We can express the KL divergence between qx and p as 



KL[ qi \\p] = E^log 



qifihH I A) 



p(fil:H\V,C, Q). 



(9) 
(10) 



where we define 



I A) 



(11) 



.p(fil:H,V\C, Q) 

Here, E gA denotes an average over using qifii.H I A). Because KL [g^ || p\ is non 



negative, (10) implies that, for all distributions qx, 



£(A;C,n)<logp(P|C,n). (12) 

The likelihood is called the "evidence" in some contexts; we therefore christen C(A; £, Q) 
the evidence lower bound (ELBO) function. 

Using ( fTO] ), we can formulate a maximization problem equivalent to ([8]), having the 
same optimal point qifii.H I A*): 



A* = argmax f, Q) . (13) 



The equivalence of @ and ([13]), together with the bound (12), shows that the best approxi 



mation in Q to the posterior yields the tightest lower bound on the marginal likelihood. 



We work with ( fT3| ) rather than ([8]) - to evaluate the KL divergence, we would need to 



compute the marginal likelihood, bringing us back to our original problem. The variational 



parameters to be adjusted are X = {h\-.h, ^v.h}- We conduct the variational inference ( [13] ) 
using block coordinate ascent on the coordinate blocks ju\, Si, juh, The coor- 
dinate updates do not have a closed form, but each update solves a smooth, unconstrained 
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convex optimization problem, as we show in Appendix [Aj There we also give a closed- 
form gradient and Hessian for the juh update, as well as closed-form gradients for the 
update under two different parametrizations. 

To finish with empirical Bayes, we explain how to obtain approximate MLEs £ and Q. 
Notice the variational inference procedure ( [T3] ) yields an optimal value A* for fixed f and C2. 
We alternate this variational inference step with a complementary optimization over £" and 
C2. In fact, these optimizations constitute a version of the expectation-maximization (EM) 
algorithm for computing MLEs. The standard E-step, where we compute the posterior 
expected complete log likelihood, is replaced with a variational E-step, where the expected 
complete log likelihood is approximated using qifti.H I X*). The variational EM algorithm 
alternates between the following steps until Q, and the variational parameters stabilize: 

E-step (variational). Using the current f and Q, run the block coordinate ascent algorithm 
as described in Appendix [Aj yielding new variational parameter values {n\. H , 

M-step. Using the current variational parameter values, update the empirical Bayes pa- 
rameter estimates: (if, Q) <— argmax^ Q C{[x\. H , C> Q). 

It can be shown that the variational E-step finds a new value of X which moves the lower- 
bounding function £(A; f , Q) towards the true log-likelihood log p(V | f, Q) from below. 
The M-step maximizes this adjusted lower-bounding function over Q and Q, as a surro- 
gate for the true log-likelihood. We then re-tighten and re-maximize until convergence. 
Appendix [A] gives details on initialization and the M-step update. 

3.2 Variational hierarchical Bayes 

Fully Bayesian inference in the mixed multinomial logit model requires calculations under 
the posterior p(fi\,H, ( , Q I V>) given in ([7]). In one sense, the setting is simpler than empir- 
ical Bayes: there are no unknown top-level parameters to estimate. All we need is to extend 
the previous section's variational inference procedure to include f and Q. Appendix |A.4 
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reports the details behind the extension; here we summarize the main ideas. 

Although the joint posterior ([7]) is not factorized, we continue to use a family Q of 
factorized distributions for the variational approximation: 



H 



Q 3 q(Pv.H, C, a I A) := q(c I H , S f ) q(n | T -1 , coj \[q{fi h \ M h, S A ) . (14) 

/!=1 

Using a factored family for a non-factored posterior is commonly called mean-field varia- 
tional inference. In ( fT4] ), | is a A^-variate normal density; q(Q | Y , co) is an 
inverse Wishart density; and the <?(/?/,) factors are ^-variate normal densities as before. In 
the analysis and the algorithm, it is convenient to use a well-known equivalence, treating 
q(Q. | Y , co) as a Wishart distribution W(Y, co) on QT . We are therefore optimizing 
over variational parameters X := [ju^, E^-, Y, co, The variational problem for 
hierarchical Bayes is to find the best approximating distribution qifii-.H, £, ^ I A*) in Q by 



solving the analog of ( 13 ): 



X* = argmax C{X) . (15) 

IsA 

As with empirical Bayes, we use a block coordinate ascent optimization algorithm to 
solve < [T5] ), iterating through the coordinate blocks that define X. Here again, all coordi- 
nate updates are convex optimizations. The details appear in Appendix [Aj updates for ju^, 
Y all have simple closed forms; co has a closed form which requires no updating; and 
the jUh and £/, updates are similar to the empirical Bayes case. 

4 Empirical results 

We compared the accuracy and speed of the variational methods described in the previous 



section to a standard and widely used MCMC approach (Allenby and Rossi 2003), on a 



suite of simulated data sets. Each data set was generated using the discrete choice model 
given by ([3]) and (Q. To simulate a data set with J choice items, K item attributes, and 
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H agents, we first fixed values of £" and Q, the parameters controlling the distribution of 
preferences in the agent population. We then independently drew a /?/, vector for each 
agent, according to @. We also drew for each agent 25 iid 7 x K item attribute matrices 
Xht consisting of iid N{0, 0.5 2 ) entries. Finally, for each agent, we used Xht and /?/, to 
simulate 25 choice events yi lt , according to ([3]). Thus, in our data sets, each agent has 25 
observed choices. 

We simulated a total of 32 different scenarios by varying 7, K, H, and the selection of 
£■ and Q. Specifically, each data set corresponds to a distinct configuration of the following 
candidate values: 3 or 12 choice items 7; 3 or 10 item attributes K; 250, 1000, 5000, or 
25000 agents H; and "low" or "high" heterogeneity of the agent population. In the low- 
heterogeneity scenario, the ^ x 1 vector £" consists of evenly spaced values from -2 to 
2, and the K x K matrix Q. is 0.25 times the identity matrix. In the high-heterogeneity 
scenario, £ is the same, but Q. is the identity matrix. The data sets with high heterogeneity 
have much more diverse collections of preference vectors fa . 

We ran variational empirical Bayes (VEB), variational hierarchical Bayes (VB), and 
the standard MCMC algorithm on the observable data from each of the 32 simulation sce- 
narios. For VEB, we declared convergence as soon as an E-step/M-step iteration caused 
the parameter estimates' joint Euclidean norm to change by less than 10 -4 , relative to their 
norm at the beginning of the iteration (here we mean the joint norm of all the variational 
parameters, together with the model parameter estimates £ and Q). The convergence crite- 
rion for VB was the same, except £ and Q. do not have point estimates - we look instead 
at the change in the variational parameters corresponding to their posterior approximation, 
namely ju^, and T. 

Choosing MCMC convergence criteria is more delicate. We tried to set the number of 



burn-in iterations and the thinning ratio algorithmically, using the technique of Raftery and 



Lewis (1992). But on several of our data sets, typical control parameter values, such as 

11 



the default settings for the raf tery . diag function in the R package coda (|Plummer 



et al. 2006), led to a very large number of burn-in iterations. Trace plots of the sampled 
parameters indicated these large burn-in values were unnecessary, so using them would 
have been unfair to MCMC in our timing comparisons. Instead, we manually investigated 
MCMC convergence and autocorrelation for several data sets, using trace plots, partial 
autocorrelation functions, and related diagnostics available in the coda package. Based 
on these studies, we chose to use 1,000 iterations of burn-in and a 10:1 thinning ratio. On 
each data set, we therefore ran 6,000 total iterations of MCMC to generate 500 draws for 
the approximate posterior. These numbers are as small as we could reasonably make them, 
in order to be fair to MCMC in the timing comparisons. 

4.1 Accuracy 

Our measure of accuracy for each inference procedure is based on the predictive choice 
distribution. Informally, this distribution gives the item choice probabilities exhibited by 
the "average" agent, when shown an item attribute matrix x new - In our simulations, we 
know the true values of Q and Q., so we can compute the true predictive choice distribution: 



new I ^newi fi)p(fi\C,a)dfi (16) 
= Epp(y new | x new ,fi) . (17) 



Equation ( |T7[ ) explains the "average agent" interpretation of the predictive choice distribu- 
tion. A slightly different take on (fT6])-([T7]) is the following: if we want to forecast the item 
choice probabilities of a new, previously unobserved decision-maker, we can think of her 
as a random draw from the agent population. Under our model, the choice probabilities for 



such a randomly drawn agent are precisely ( 16 )-( 17 ) 



For any particular x new , we use the true values of £" and Q. to compute a Monte Carlo 
estimate of ( fTTj ). We base the estimate on enough draws of /? ~ M(£, O.) to insure that 
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its variability does not affect even the least significant digit of our reported results. We 
handle the integral over /? in the same way for the estimated predictive choice distributions 



furnished by each of the three inference procedures. For VEB, the estimate is ( 16), with f 



and Q. replaced by £veb an d £&veb- O n me other hand, with VB and MCMC, we obtain a 



posterior distribution over £" and Q; we take the mean of ( |16| ) under this posterior as a point 
estimate of the predictive choice distribution: 



iHVnew I -''•new > 



Piy^sw I -''•new 5 



p(C,Q\V)dCdQ . (18) 



For VB, the posterior density p(£, Q. \ V) in ( [18] ) is approximated by the fitted variational 
distribution 

q(C\fi l: ,^)q(a\T-\co) . (19) 
For MCMC, the posterior is approximated as usual by the empirical distribution of draws 



from a Markov chain. In both cases, we handle the integral over f and Q. in ( fl8| ) with 
another exhaustive Monte Carlo approximation. 

We measure the error of each inference procedure as the distance from its estimate of the 
predictive choice distribution to the true distribution. As the metric on distributions, we use 
total variation (TV) distance, leading to what we call the "TV error" of each procedure. In 
this setting, TV error equals the maximum, over all choice-item subsets, of the difference in 
the probabilities assigned to the subset by the estimated versus the true choice distribution. 
We also need to choose the attribute matrix jc new at which the true and estimated predictive 
choice distributions are calculated. In each simulation scenario, we computed the TV error 
of VEB on 25 different random draws of x new . We then compared the three procedures 
using the x new which yielded the median of the 25 TV errors. In this sense our results are 
representative of accuracy for a "typical" item attribute matrix. However, results using any 
one of the 25 matrices were qualitatively the same in the examples we checked. 
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3 attributes 10 attributes 



Agents 


Het. 


VEB 


VB 


MCMC 


VEB 


VB 


MCMC 


250 


Low 
High 


0.23 
0.49 


0.31 
0.48 


0.35 
0.92 


0.66 
0.68 


0.50 
0.51 


0.59 
1.46 


1,000 


Low 
High 


0.74 
0.53 


0.63 
0.33 


0.46 
0.09 


0.49 
0.30 


0.17 
0.39 


0.14 
0.18 


5,000 


Low 
High 


0.25 
0.19 


0.37 
0.25 


0.35 
0.39 


0.56 
0.35 


0.59 
0.20 


0.07 
0.38 


25,000 


Low 
High 


0.63 
0.53 


0.74 
0.53 


NA 
NA 


0.59 
1.60 


0.55 
1.10 


NA 
NA 



Table 1 : Total variation error in percentage points, for simulated data sets with three choice 
items. MCMC results are unavailable in the 25,000 agent case because the sampler ex- 
hausted memory resources before converging. See the text for the definition of total varia- 
tion error. 



TV error values for the three inference procedures are presented in Table [T] for the 3- 
item simulation scenarios, and in Table [2] for the 12-item scenarios. The main conclusion 
we draw from Tables [T]and|2]is simple: on these data sets, there are no practical differences 
in accuracy among VEB, VB, and MCMC. The scale of the TV error for all the procedures 
is the same; that scale is larger in the 12-item case than the 3-item case, but all three 
procedures exhibit high accuracy on all data sets. The magnitude of our simulation study 
makes it difficult to carry out the replications required to put standard errors in these tables. 
But even if the differences among TV errors in every scenario were "significant" under a 
suitable definition, the patternless alternation in the identity of the most accurate method 
would make more specific conclusions dubious. 

Figure [T] shows in a different way the comparable accuracy of the three procedures. 
When there are three choice items, any predictive choice distribution can be plotted as a 
point in the triangular simplex. The figure shows the close proximity of each procedure's 
estimated choice distribution to the true distribution in one simulation scenario. Plots for 
all the other three-item scenarios are qualitatively the same. Figure [T] also shows contours 
of the VB and MCMC approximate posterior distributions. We see that VB is producing 
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3 attributes 



10 attributes 



Agents 


Het. 


VEB 


VB 


MCMC 


VEB 


VB 


MCMC 


250 


Low 
High 


o o o 
z.oo 

1.44 


O OA 

2.80 
1.94 


O £ A 

2.o4 
1.64 


1.97 
2.43 


1.91 

2.37 


O A A 

2.44 
2.62 


1,000 


Low 
High 


1.11 
0.98 


1.27 
1.18 


1.09 
1.18 


0.99 
1.99 


1.00 
2.05 


1.60 
1.96 


5,000 


Low 
High 


1.25 
1.14 


1.45 
0.98 


1.18 
1.10 


0.95 
0.71 


1.13 
0.91 


0.97 
0.92 


25,000 


Low 
High 


0.22 
0.99 


0.33 
0.57 


NA 
NA 


0.51 
1.23 


0.53 
0.96 


NA 
NA 



Table 2: Total variation error in percentage points, for simulated data sets with 12 choice 
items. See the caption accompanying Table [T] 




Figure 1: Triangle plot of the true predictive choice distribution and its estimates in the 
three-item case, for the simulation with 250 decision-makers, 10 attributes, and high het- 
erogeneity. See accompanying text. 



15 



not only a posterior mean similar to MCMC, but also a similar posterior density in the 
neighborhood of the mean. 

4.2 Speed 

For each simulated data set, we ran the three procedures in turn on the same unloaded 
machine: a 64-bit dual-core 3.2 GHz Intel Xeon processor with 8 GB of main memory. 
For the MCMC inference, we used the rhierMnlRwMixture function in the R package 
bayesm ( |Rossi and McCulloch||2007"| ), which has efficient vectorized implementations of 
the inner sampling routines. This package stores all MCMC draws in memory, however. 
For our largest data sets, with 25,000 decision-makers, the machine's memory was ex- 
hausted before MCMC converged. We were able to run MCMC for 1,000 iterations in this 
case, which allowed us to extrapolate accurately the time that would have been required 
for 6,000 iterations. We implemented the variational algorithms in R, with compiled C 
components for the numerical optimization routines. 

Figure [2] displays time to convergence on each data set for the three procedures, accord- 
ing to the convergence criteria previously described. Within each panel, convergence time 
is plotted as a function of the number of agents, for fixed values of the other simulation 
parameters. Note that the vertical axis shows convergence time on a logarithmic scale, to 
ease comparison of MCMC to the variational methods. All the procedures scale roughly 
linearly with the number of agents, which leads to the logarithmic curves seen in the figure. 
The conclusions are the same in all the scenarios we simulated: variational methods con- 
verge faster than MCMC, and the magnitude of the difference increases with the number 
of observed agents. MCMC uses two days of computation time for 25,000 agents with 
12 choice items, 10 item attributes, and high heterogeneity, versus an hour for each of the 
variational techniques. In the same setting but with low heterogeneity, MCMC's two-day 
computation compares with two hours for VEB and six hours for VB. In other scenarios, 
variational run times are measured in minutes, as opposed to hours or days for MCMC. 
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Figure 2: Timing results for variational empirical Bayes (VEB), variational hierarchical 
Bayes (VB), and MCMC. Within each panel, convergence time is plotted on the log scale 
as a function of the number of agents, for fixed values of the other simulation parameters 
(shown at the top of each panel). 
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5 Discussion 



Variational methods allow estimation of hierarchical discrete choice models in a small 
fraction of the time required for MCMC. They open Bayesian inference to a wide class 
of problems for which resource constraints make MCMC intractable, such as the MML 
model with many heterogeneous units. For now, variational methods appear to be the only 
viable option in these cases. 

Of course, one can use variational methods to estimates many more types of models 
than the MML model examined here. Within the MML family, it would be straightforward 
to add a utility scaling parameter, or allow the heterogeneous coefficients themselves to 
depend on observed covariates. The value of the variational approach is greatest when sub- 
sampling of data is ill-advised. For example, consider a linear model with heterogeneous 
coefficients on exogenous and endogenous covariates, where the available instrumental 
variables only weakly explain the endogenous part. To draw inferences about the covari- 
ances of the heterogeneous parameters, we may need a large amount of data to achieve 
reasonable power in hypothesis testing. MCMC is untenable here, but variational methods 
have promise. When factorized variational distributions are inadequate, alternatives such 
as mixtures of normals or Dirichlet processes ( Blei and Jordan|20 06) can be applied. 



We emphasize that we do not advocate abandoning MCMC in favor of variational meth- 
ods. On the contrary, we suggest using MCMC when possible. MCMC offers consistency 
guarantees with no analog in variational methods, and it has withstood decades of scrutiny. 
But we advise against subsampling the data (i.e., throwing out information), or discarding 
key modeling elements, simply to make the problem fit within the time and resource con- 
straints of MCMC. The possibility of applying variational methods to previously intractable 
problems makes them an important addition to the statistician's toolkit. 
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A Variational inference and parameter estimation 

In this appendix we describe the variational inference and estimation procedures for the 
mixed multinomial logit model. A few words on notation: A y means the matrix 
A is positive definite; A y means positive semidefinite; |A| is the determinant of A. 
For a scalar function s : K — > K and a vector v e K n , s(v) means the n x 1 vector 
OOi), . . . ,s(v n )) J . 

A.l The empirical Bayes ELBO 



For empirical Bayes in the MML model, the ELBO objective function (11) becomes 



H H T h 

H(q) + X E ? lo § P(fih I f , + X X E ? lo § 1 • (20) 

h=\ h=l t=l 



The first term in ( |20| ) is the Shannon entropy of the variational distribution. The second and 
third terms are (minus) an unnormalized cross entropy - the missing normalization constant 
is the marginal likelihood. Recall that the variational distribution qifti.H I Mi.H, ^v.h) is a 
product of normal distributions Af(jUh, 



The first and second terms of pO] ) are straightforward to derive; the third term requires 
more attention. Using the multinomial logit mass function 



pijht I Xkt,flh) = \ } 



(21) 



the third term becomes 



H T h 

ZZ 

h=l t-\ 



^ylMhtj^h)-^ log 



(22) 



The expected log-sum-exp in ( |22| ) has no closed form. For variational inference, we there- 
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fore approximate the ELBO objective function C using a new objective function C. To 
construct C, we consider two alternatives: the zeroth-order and first-order delta method for 
moments ( Bickel and Doksum|2007 ), which we call DO and Dl respectively. DO is equiva- 
lent to applying Jensen's inequality to the expected log-sum-exp, resulting in the following 



lower bound to (22): 



H T h 



[DO] 



h=l t=l 



^y J ht( x htjMh) ~ log( ^expO^,, + (\/2)xl tj ^ h x htj )\ 



. (23) 



Here we used the usual formula for the mean of a lognormal random variable. In a different 
context, 



Blei and Lafferty (2007) consider an approximation equivalent to DO but expressed 



using a redundant variational parameter. 



For approximation Dl, we restrict Z/, y to be diagonal, and define 
ah := log(diag{£/J) e M. K . Using results in Appendix[Bj we obtain the following approx- 
imation to ([22]): 

^y J ht( x htj{*h) -lognTexpCx^/Oj - ^exp(<7/,) T (///,) 

(24) 



™ EE 

/7=1 /=! 



with &(jUh) e as defined in Appendix |b} Notice that, unlike DO, approximation Dl 
does not preserve the guarantee that the optimal value of the variational optimization lower 
bounds the marginal likelihood. However, in our simulations, using Dl resulted in more 
accurate variational approximations to the posterior. 



In this appendix we give a derivation based on approximation DO. The derivation for 
Dl is similar, but simpler, because Z/, is treated as diagonal. Under DO, the final empirical 
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Bayes objective function is 



C(M hH, C, = - jr log [(27T g) z I *h I] 

- I log ((2t0* |Q|) - i tr jo" 1 £ { I, + - - cf} J 



(25) 



+ZZ 



The first line in ( [25] ) uses the well-known entropy of the normal distribution. The second 
line uses the cross-entropy of two normal distributions, also well known. The third line is 
approximation DO. 

A.2 Empirical Bayes variational E-step 



Here we describe a block coordinate ascent algorithm to maximize ([25]) over the varia- 
tional parameters ju\-h and Although the problem is not jointly convex in all these 
parameters, each fih and E/, coordinate update solves a smooth, unconstrained convex opti- 
mization problem. The requirement S/, >z is satisfied after each update. We initialize the 
variational parameters at the maximum likelihood estimates from a homogeneous model 
(in which all agents share a common value). 



The concavity of ( [25] ) in juj t follows from the fact that Q. y and from the convexity 
of the log-sum-exp function. We update fih using standard algorithms for unconstrained 



convex optimization ( |Boyd and Vandenberghe| |2004), supplying an analytic gradient and 
Hessian as follows. Define the function w(ju, E, x) taking values in ~R. J , with jth compo- 
nent exp (xj [x + (1/2)jcJSjC;V and normalized to sum to one across j. The gradient of 
C with respect to jXh can then be written 



dfx h 



T h J 



= -q i (ju h - c) + y, y ^yj t - w j (ju h , Zh, xht)^xht 



(26) 



t=i j=i 
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Note the similarity of this gradient to the gradient from an L 2 -regularized multiple logistic 
regression: it consists of a contribution from the regularizer (the left-hand term), plus a 
residual-weighted sum of covariate vectors. Abbreviating w(ju.h, Xht) to wj 7t , an argu- 
ment using matrix differentials (Magnus and Neudecker|2007 ) gives the Hessian 



dfihdju 



Y = -&■ 1 - [*I dia 8{%}*to _ ( x ht w ht)( x ht w ht) T ] • (27) 



The E^ coordinate update is harder, because we need to insure that E& >z 0. Using a 
reformulation, we can avoid making the constraint explicit, which would complicate the 
optimization. Let E/, = LhLj for a lower-triangular matrix L^. Since E/, >z 0, one such 
Lf, always exists — the Cholesky factor. We replace each E/; in C with E/^L/,) '■= LhLj, 
and optimize over the unconstrained set of lower- triangular matrices Lf,. 

The objective function ( |25] ) remains concave in L/,. To see this, compare the terms 
depending on £/, = LhL^ to the function studied in Appendix|cj We now give the gradient 
with respect to L/,. Standard matrix differentiation of p5| ) leads to the E/, gradient 



5£ 1 



92* 2 



7ft 



(28) 



Again using matrix differentials and the Cauchy invariance rule, it is not hard to show that 
the gradient with respect to is 



dC 
dL h 



'■(j^j L h = L h T ~ + Y^xldmg{w ht }x h ^jL h . (29) 



Note that this is the gradient with respect to a dense matrix Lh- Since we optimize over 
lower-triangular matrices, i.e. vech(L^), we need only use the lower triangular of the gra- 
dient. This is convenient for the term L^ 1 : it is upper-triangular, so its lower triangle is a 
diagonal matrix. Furthermore, from a standard result of linear algebra, the diagonal entries 
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are simply 1 ft a, where the £a's form the diagonal of L^. 

In practice we do the and E/, updates in a single step by optimizing jointly over juh 
and Lh, which remains a convex problem. 

A.3 Empirical Bayes M-step 



In the M-step, we maximize ( [25] ) over Q and Q. Identifying the terms which depend on £, 
we recognize the usual Gaussian mean estimation problem. Further, ((25]) is easily seen to 
be concave in Q.~ l , with a closed-form solution of the corresponding first-order condition. 
We obtain the M-step updates 

I H 1 H _ 

C<-Jj^Mh, Q <- + (30) 

h=\ h=l 

Here Cov(//.) is the empirical covariance of the fi^ vectors. 
A.4 Variational hierarchical Bayes 

In the fully Bayesian MML model, £" and Q. have prior distributions, with corresponding 
variational factors given in ( ]T4| ). The ELBO in this case has the same form as ( f20j ), with 
two differences. First, H(q) contains two new terms 

H (q(c | H , + H(q(Q \ T _1 , to)) . (31) 

Second, there are two new cross-entropy terms 

log p(C I y^o, Q ) + E 9 log p(Q | 5, v) . (32) 



Also, the middle term of ( |20| ) changes in the fully Bayesian case, because Q and Q are now 
averaged over rather than treated as constants. 

Using known formulas for normal and Wishart entropies, the two new entropy terms 
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are seen to equal 



1 T , v co — K — 1 coK , \ 

- log [{2ne) K \Y. c \] D(a>, T) + — + A ra (T) 



(33) 



Here we used the expected log determinant of a Wishart random matrix 



D(co, T) := log (2 K \T\) + ]T> ( ~ + ^ - ) 



(34) 



and the log normalization constant of the Wishart distribution 



A ra (T) := log 



2<oKp. K(K-\)/A 



E[r(^)' 



CO 

-K — log |TT| 



(35) 



(see, for example, [Beal|2003] ). The new cross entropy terms for f and Q work out to 



- l - {log [(27r)*|n |] + tr (Qq 1 [S f + (fi ( - fi )(ji c - /? ) T ]) } (36) 



and 



, v — — 1 , CO / i v 

- A^S- 1 ) + D(a>, T) - - tr^T) 



(37) 



respectively. The middle term of ( 20 ) eventually becomes 



-y{Zlog(2*)-D(G>,T)} - 



£0 



tr 



\t(hX c + J^(U + (Mc-Mh)(M{-Mh) T )^ ■ (38) 



With these changes, it is not hard to see that £ is concave separately in ju^ and £ 4 -. The 
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first-order conditions for block coordinate ascent lead to the updates 



H 



H <- (a^ + HcoT) MQ-Vo + ^T^^V (39) 

^ h=i ' 

Z f <- (Q.- 1 + HcoT)~ l . (40) 

By inspection, h 0, so this constraint need not be explicitly enforced. Note the similar- 
ity to conjugate posterior updating: on the precision scale, E f is the sum of the prior pre- 
cision matrix 1 and H copies of the variational posterior mean coT for Similarly, 
fi£ is a precision- weighted convex combination of the prior vector /?o and the empirical 
average of the variational posterior means fii-.H for 

The updates for T and co are similarly straightforward to derive; we obtain 

co 4 — v + H , (41) 
T<- ^S-'+^Xh + iMz-MhXMz-jUh^+HX^ . (42) 



Notice that the solution (41 ) for co involves only the constants v and H . We compute co 



once in advance, leaving it unchanged during the variational optimization. 

B An application of the delta method 

Let f(v) be a function from H K to R. According to the multivariate delta method for 
moments (Bicke l~a"nd Doksum|2007[ ), 



Ef(V) * f(EV) + \ tr [(f^r) Cov(V) 



(43) 



Consider the case 

/(o)=log(i T exp(xo)) , (44) 
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where x is a J x K matrix whose rows are the vectors xj. Let V ~ J\ffc(pi, £)> an d restrict 
Z such that E = diag{exp(cr)} for a e M. K . We can now rewrite d43l): 



Elog(l T exp(xV)) w log(l T exp(jc/i)) + ^0(/^) T exp((r) , (45) 



where 0(//) is the diagonal of the Hessian of /, evaluated at the point pi. Define s = 
1 T exp(x/u). Using matrix differentials, it can be shown that 

&(pi) = s~ l (x O x) T exp(xpi) — s~ 2 (x T exp (*//)) (x T exp(xpi)) , (46) 

where denotes the Hadamard product. 

To use the approximation ( [45] ) in an optimization over pi, we need to compute the 
gradient. The formula for makes this a more extensive but still mechanical exercise 
in differentials. One obtains 



5 _ 1 T xii T 

— = s l x 1 e Xfl + -x 1 
dpi 2 



(s~ l diag {e*P} - s~ 2 e x » (e x *) J \ (x O x) + 



2(s~ 3 e Xfl {(x T e XfI ) © (x T e Xfl )} T - s ~ 2 diag {e x ' i } x diag {xV^})] exp(a) . (47) 

C A convexity result 

Let a\, . . . , ad be scalars, ci, . . . , Cd be n-vectors, p, r > 0, and 2^0. We show here 
that the function 



f(B) = r\og\BB T \ -ptr(QBB T ) - log/^exp [aj + c]BB t cA 



(48) 



is concave on the set of full-rank n x n matrices. 

We argue that each of the three constituent terms, from left to right, is concave. The 
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second differential of g(B) = r log! BB T I is 



d 2 j? = dtr{2rfl _1 dfl} = tr 



-2r 



B 



(49) 



By Theorem 10.6.1 of (Mag nus and Neudecker|2007l ), the Hessian of g is 

—2rK n (B~ T (g> where K n is the order-n commutation matrix and (g> denotes the 

Kronecker product. We now show that K n (B~ T ® is (matrix) positive-definite. 



(dvecX) 7 ^ (V T <g> (dvecX) = (dvecX) 7 ^ vec {B _1 (dX)5 -1 } 

= (vecdX) T vec{5- T (dX) T 5- T } 
^5 _1 dx) 2 J > 0. 



= tr 



(50) 
(51) 
(52) 



Equation ( [50| ) follows from the well-known fact that vec ABC = (C T <g> A) vec B. Thus, 
the Hessian of g is negative definite, and r logl BB T I is concave. 



Concavity of the middle term in ( [48] ) follows in the usual way from the univariate con- 
vexity of the function 



g(t) := tr (Q(M + tP)(M + tP) T ) = ^(m ; + t Pi ) T Qim, + t Pi ) (53) 



for fixed matrices M and P, with columns m ; and p, . To see that the rightmost term in ( [48] ) 
is concave, define 

gj (t) := aj + c](M + tQ)(M + tQ) T Cj 
for j = 1, . . . , d and fixed matrices M and Q. Each g ; is convex in t, and the rightmost 



term in ( [48] ) is (minus) the log-sum-exp function composed with the gj's. Concavity of this 
term in t, and hence in B, follows from ( |Boyd and Vandenbe rghe|2004] p. 86). 
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